

#include "glodef.h"
#include "errlog.h"
#include "f3io.h"
#include "vecs.h"
#include "glovar.h"
#include "grid3d.h"

#include "soah_common.h"

#include <stdio.h>
#include <stdlib.h>
#include <string.h>


int MAX_SHAPES = 10000;

static char saved_log[FILENAME_MAX];
static char saved_root[FILENAME_MAX];
static char saved_log_dir[FILENAME_MAX];
static int shape_count;
static int DONE;
static node_t *shape_nodes;
static vector init_u_c;
static scalar init_p_c;
static vector finl_v_c;
static scalar finl_q_c;
static double sg_norm_init;
static double sgn_c, sgni_c;
static scalar sg_xi_c;

static void soah_shape_optimize_begin(vector init_u, scalar init_p)
{
    soah_check_heart_geometry();
    alloc_sg();
    make_sg_z();
    clean_scal(sg_xi);

    GetLogfilePath(saved_log, FILENAME_MAX);
    strncpy(saved_root, root_dir(), FILENAME_MAX-1);
    strncpy(saved_log_dir, logs_dir(), FILENAME_MAX-1);

    init_u_c=alloc_vec(get_N_vec(init_u));
    copy_vec(init_u, init_u_c);
    init_p_c=alloc_scal(get_N_scal(init_p));
    copy_scal(init_p, init_p_c);
    finl_v_c=alloc_vec(get_N_vec(init_u));
    finl_q_c=alloc_scal(get_N_scal(init_p));
    sg_xi_c=alloc_scal(get_N_scal(sg_xi));

    if((shape_nodes=malloc(grid->NN*sizeof(node_t)))==NULL)
        error_msg("%s: memory allocation error.\n", __func__);
    memcpy(shape_nodes, reference_nodes, grid->NN*sizeof(node_t));

    // init some global parameters
    shape_count=1;
    sg_lambda=0.0;
    sg_norm=0.0;
    sg_norm_inf=0.0;
    sg_norm_init= -1.0;
    DONE=FALSE;
}

static void soah_shape_optimize_end(void)
{
    free_sg();
    free_vec(init_u_c);
    free_scal(init_p_c);
    free_vec(finl_v_c);
    free_scal(finl_q_c);
    free_scal(sg_xi_c);
    free(shape_nodes);
}

/****************************************************************************/

/****************************************************************************/

static void start_shape_log(int shape_count, char *shape_log)
{
    static char shape_dir[FILENAME_MAX];
    sprintf(shape_dir, "%s/shape_%i", saved_root, shape_count);
    set_root_dir(shape_dir);
    create_dir(root_dir());
    create_dir(logs_dir());
    sprintf(shape_log, "%s/shape_%i.log", saved_log_dir, shape_count);
    PrintLog("Logging to file %s\n", shape_log);
    InitLogfile(shape_log);
    open_dir(root_dir());
}

/****************************************************************************/

/****************************************************************************/

static int create_shape_mesh(const char *shape_log)
{
    static char mesh_log[FILENAME_MAX];
    int ret;
    sprintf(mesh_log, "%s/mesh.log", logs_dir());
    PrintLog("Logging mesh generation to file %s\n", mesh_log);
    InitLogfile(mesh_log);
    if((ret=soah_move_solid_wall_mesq())==FALSE)
        return FALSE;
    ret=soah_generate_meshes(FALSE);
    AppendLogfile(shape_log);
    if(!flags.moving_grid)
    {
        if(massP) fe_matr_assemble(grid, 1, &massP);
        if(massU) fe_matr_assemble(grid, 1, &massU);
    }
    return ret;
}

/****************************************************************************/

double soah_shape_cost(vector init_u, scalar init_p)
{
    int ret;
    switch(flags.periodic)
    {
        case per_no:
            ret=solve_state_ivp(init_u, init_p);
            break;
        case per_march:
            ret=solve_state_periodic(init_u, init_p);
            break;
        case per_adjoint:
            ret=solve_state_periodic_adjoint(init_u, init_p);
            break;
        case per_armijo:
            ret=solve_state_periodic_armijo(init_u, init_p);
            break;
        case per_broyden:
            ret=solve_state_periodic_broyden(init_u, init_p);
            break;
        default:
            error_msg("%s: flags.periodic unrecognized.\n", __func__);
    }
    if(!ret)
    {
        total_cost=1e99;
        warning_msg("Failed to solve state problem.  "
            "Setting cost to %f to force step rejection.", total_cost);
    }
    return total_cost;
}

static int soah_shape_gradient(vector finl_v, scalar finl_q)
{
    int ret;
    switch(flags.periodic)
    {
        case per_no:
            // clean_vec(finl_v);
            // clean_scal(finl_q);
            ret=solve_adjoint_shape_fvp(finl_v, finl_q);
            break;
        case per_adjoint:
        case per_armijo:
        case per_lbfgs:
            PrintLog("%s: this flags.periodic value is not available for the "
                "adjoint problem.\n", __func__);
        case per_march:
            ret=solve_shape_adjoint_periodic(finl_v, finl_q);
            break;
        case per_broyden:
            ret=solve_adjoint_shape_periodic_broyden(finl_v, finl_q);
            break;
        default:
            error_msg("%s: flags.periodic unrecognized.\n", __func__);
    }

    if(!ret)
        warning_msg("%s: failed to solve the adjoint problem\n", __func__);
    else
        compute_shape_gradient();

    return ret;
}

/****************************************************************************/

static int soah_shape_step_begin(void)
{
    char shape_log[FILENAME_MAX];

    PrintLogLine(78, '=');
    PrintLog("Attempting shape with sg_norm = %g and sg_lambda = %g\n",
        sg_norm, sg_lambda);

    PRINTF("Attempt shape with lambda = %g\n", sg_lambda);
    start_shape_log(shape_count, shape_log);

    if(!create_shape_mesh(shape_log))
    {
        AppendLogfile(saved_log);
        PrintLog("Failed to create mesh.\n");
        return FALSE;
    }
    save_ref_grid();
    PrintLog("GridVolume is %.8g\n", GridVolume(grid));
    return TRUE;
}

static void soah_shape_step_accept(vector init_u, scalar init_p,
    vector finl_v, scalar finl_q)
{
    copy_vec(init_u, init_u_c);
    copy_scal(init_p, init_p_c);
    copy_vec(finl_v, finl_v_c);
    copy_scal(finl_q, finl_q_c);
    copy_scal(sg_xi, sg_xi_c);
    sgn_c=sg_norm;
    sgni_c=sg_norm_inf;

    AppendLogfile(saved_log);
    PrintLog("       %10s   %-20s = %.8g\n", "total_cost", "", total_cost);
    PrintLog("       %10s  %-20s = %.8g\n", "sg_norm_inf", "", sg_norm_inf);
    PrintLog("       %10s   %-20s = %.8g\n", "shape_grad", "", sg_norm);
    PrintLog("Passed acceptance check.  Accepting step.\n");

    PRINTF("\rshape %u : total_cost=%7g sg_norm=%7g sg_norm_inf=%7g\n",
        shape_count, total_cost, sg_norm, sg_norm_inf);

    // suggested initial lambda for next shape step
    sg_lambda=grid->MAXH*0.3/sg_norm_inf;
    memcpy(shape_nodes, reference_nodes, grid->NN*sizeof(node_t));
    shape_count++;
}

static int soah_shape_step_acceptance_check(double F, double Fc)
{
    const double alpha=1e-3;
    // acceptance check
    if(F-Fc > -alpha*sg_lambda*sg_norm_inf)
        return FALSE;
    else
        return TRUE;
}

static int soah_shape_step_termination_check(void)
{
    if(shape_count > MAX_SHAPES)
        return (DONE=TRUE);
    const double ta=1e-6;
    const double tr=ta*1e3;
    if(sg_norm_init<=0.0) sg_norm_init=sg_norm;
    if((sg_norm<(tr*sg_norm_init+ta))&&(shape_count>3))
        return(DONE=TRUE);
    else
        return(DONE=FALSE);
}

/****************************************************************************/

/****************************************************************************/

/****************************************************************************/

void soah_shape_optimize(vector init_u, scalar init_p,
    vector finl_v, scalar finl_q)
{
    int DONE=FALSE;

    soah_shape_optimize_begin(init_u, init_p);

    shape_count=0;
    while(!DONE)
    {
        static char shape_log[FILENAME_MAX];

        shape_count++;

        start_shape_log(shape_count, shape_log);

        create_shape_mesh(shape_log);

        /* initial conditions */
        clean_vec(init_u);
        clean_scal(init_p);

        soah_shape_cost(init_u, init_p);

        clean_vec(finl_v);
        clean_scal(finl_q);
        error_msg("%s: this function does not work.\n", __func__);
        // solve_shape_adjoint_periodic(finl_v, finl_q);
        // solve_adjoint_shape_periodic_broyden(finl_v, finl_q);

        AppendLogfile(saved_log);
        DONE=compute_shape_gradient();

    }

    soah_shape_optimize_end();

}

/****************************************************************************/

/****************************************************************************/

void soah_shape_optimize_armijo(vector init_u, scalar init_p,
    vector finl_v, scalar finl_q)
{

    const double beta=0.8;
    double Fc=1e99, F;
    double Gc=1e99, G;
    double prev_lambda=0.0;

    soah_shape_optimize_begin(init_u, init_p);

    while(!DONE)
    {
        // opens .log file and generates new mesh
        if(!soah_shape_step_begin())
            goto reject_step; // failed to create mesh

        // compute cost functional at current shape (mesh)
        F=soah_shape_cost(init_u, init_p);

        // acceptance check
        if(!soah_shape_step_acceptance_check(F, Fc))
        {
            AppendLogfile(saved_log);
            PrintLog("       %10s   %-20s = %.8g\n", "total_cost", "", total_cost);
            PrintLog("Failed acceptance check.  Rejecting step.\n");
            goto reject_step;
        }

        soah_shape_gradient(finl_v, finl_q);
        G=sg_norm_inf;
        if(G>=Gc*1.5)
        {
            AppendLogfile(saved_log);
            PrintLog("       %10s   %-20s = %.8g\n", "total_cost", "", total_cost);
            PrintLog("       %10s  %-20s = %.8g\n", "sg_norm_inf", "", sg_norm_inf);
            PrintLog("       %10s   %-20s = %.8g\n", "shape_grad", "", sg_norm);
            PrintLog("Failed acceptance check.  Rejecting step.\n");
            goto reject_step;
        }

        /* accept step */

        /* save_lambda overrides the default initial lambda
         * suggesed in soah_shape_step_termination_check */
        double save_lambda=sg_lambda;

        soah_shape_step_accept(init_u, init_p, finl_v, finl_q);
        soah_shape_step_termination_check();

        if(save_lambda>1e-9)
            sg_lambda=save_lambda;
        prev_lambda=sg_lambda;
        Fc=F;
        Gc=G;

        continue;

        /* reject step */
reject_step:
        set_reference_nodes(shape_nodes);

        sg_lambda*=beta;
        if(sg_lambda<0.001*prev_lambda)
        {
            PRINTF("lambda = %g is too small.\n", sg_lambda);
            PrintLog("lambda = %g is too small.\n", sg_lambda);

            warning_msg("Shape %d does not decrease the value of the "
                "cost functional", shape_count);
            DONE=TRUE;
            continue;
        }
        copy_vec(init_u_c, init_u);
        copy_scal(init_p_c, init_p);
        copy_vec(finl_v_c, finl_v);
        copy_scal(finl_q_c, finl_q);
        copy_scal(sg_xi_c, sg_xi);
        sg_norm=sgn_c;
        sg_norm_inf=sgni_c;
    }

    soah_shape_optimize_end();
}

void soah_shape_optimize_pit(vector init_u, scalar init_p,
    vector finl_v, scalar finl_q)
{
    double Fs;

    soah_shape_optimize_begin(init_u, init_p);

    if(!soah_shape_step_begin())
        error_msg("%s: this cannot happen.\n", __func__);
    Fs=soah_shape_cost(init_u, init_p);

    soah_shape_gradient(finl_v, finl_q);
    soah_shape_step_accept(init_u, init_p, finl_v, finl_q);
    soah_shape_step_termination_check();

    while(!DONE)
    {
        /*
         * find lmax and lmin
         */
        double lmin=sg_lambda*1e-3;
        double lmax= -1.0;
        while(soah_shape_step_begin()) // increase lambda until mesh fails
        {
            set_reference_nodes(shape_nodes);
            // set_ref_nodes(alegrid, shape_nodes);
            AppendLogfile(saved_log);
            lmax=sg_lambda;
            sg_lambda*=1.25;
        }
        if(lmax<=0)
        {
            // decrease lambda until mesh succeeds
            set_reference_nodes(shape_nodes);
            // set_ref_nodes(alegrid, shape_nodes);
            AppendLogfile(saved_log);
            sg_lambda*=0.8;
            while(!soah_shape_step_begin())
            {
                set_reference_nodes(shape_nodes);
                // set_ref_nodes(alegrid, shape_nodes);
                AppendLogfile(saved_log);
                sg_lambda*=0.8;
            }
            lmax=sg_lambda;
        }
        set_reference_nodes(shape_nodes);
        // set_ref_nodes(alegrid, shape_nodes);
        AppendLogfile(saved_log);

        PrintLog("    pit:  lmin = %g    lmax = %g.\n", lmin, lmax);
        PRINTF("  pit:  lmin = %g    lmax = %g.\n", lmin, lmax);

        double l2=sg_lambda=0.9*lmax;
        if(!soah_shape_step_begin())
            error_msg("%s: %4d: this should not happen.\n", __func__, __LINE__);
        copy_vec(init_u_c, init_u);
        copy_scal(init_p_c, init_p);
        double F2=soah_shape_cost(init_u, init_p);
        set_reference_nodes(shape_nodes);
        // set_ref_nodes(alegrid, shape_nodes);
        AppendLogfile(saved_log);
        PrintLog("    pit:  l2 = %g    F2 = %g.\n", l2, F2);

        double l0=sg_lambda=0.3*lmax;
        if(!soah_shape_step_begin())
            error_msg("%s: %4d: this should not happen.\n", __func__, __LINE__);
        copy_vec(init_u_c, init_u);
        copy_scal(init_p_c, init_p);
        double F0=soah_shape_cost(init_u, init_p);
        set_reference_nodes(shape_nodes);
        // set_ref_nodes(alegrid, shape_nodes);
        AppendLogfile(saved_log);
        PrintLog("    pit:  l0 = %g    F0 = %g.\n", l0, F0);

        double l1=sg_lambda=0.6*lmax;
        if(!soah_shape_step_begin())
            error_msg("%s: %4d: this should not happen.\n", __func__, __LINE__);
        copy_vec(init_u_c, init_u);
        copy_scal(init_p_c, init_p);
        double F1=soah_shape_cost(init_u, init_p);

        int pit_found=FALSE;
        int lmin_reached=FALSE;
        int lmax_reached=FALSE;

        while(1)
        {
            if((F0>=F1)&&(F1<=F2))
            {
                pit_found=TRUE;
                break;
            }
            if((F0<F1)&&(F1>F2))
            {
                // we are on a bump.  Which way to go?
                if(F0<F2) goto go_left;
                else /* F0 > F2 */ goto go_right;
            }
            if((F0<F1)&&(F1<F2)) goto go_left;
            if((F0>F1)&&(F1>F2)) goto go_right;
go_left:
            {
                if(l0<=lmin)
                {
                    lmin_reached=TRUE;
                    break;
                }
                set_reference_nodes(shape_nodes);
                // set_ref_nodes(alegrid, shape_nodes);
                AppendLogfile(saved_log);
                PrintLog("    pit:  l0 = %g   l1 = %g   l2 = %g\n", l0, l1, l2);
                PrintLog("    pit : F0 = %g   F1 = %g   F2 = %g\n", F0, F1, F2);
                PrintLog("    pit : going left.\n");
                PRINTF("\r  pit:  l0 = %5g   l1 = %5g   l2 = %5g  \t\t\n\t"
                    "F0 = %5g   F1 = %5g   F2 = %5g   going left\n",
                    l0, l1, l2, F0, F1, F2);
                l2=l1;
                l1=l0;
                F2=F1;
                F1=F0;
                l0=sg_lambda=0.8*l0;
                if(!soah_shape_step_begin())
                    error_msg("%s: %4d: this should not happen.\n", __func__, __LINE__);
                copy_vec(init_u_c, init_u);
                copy_scal(init_p_c, init_p);
                F0=soah_shape_cost(init_u, init_p);
                continue;
            }
go_right:
            {
                if(l2-lmax>= -1e-5*lmax)
                {
                    lmax_reached=TRUE;
                    break;
                }
                set_reference_nodes(shape_nodes);
                // set_ref_nodes(alegrid, shape_nodes);
                AppendLogfile(saved_log);
                PrintLog("    pit:  l0 = %g   l1 = %g   l2 = %g\n", l0, l1, l2);
                PrintLog("    pit : F0 = %g   F1 = %g   F2 = %g\n", F0, F1, F2);
                PrintLog("    pit : going right.\n");
                PRINTF("  pit:  l0 = %5g   l1 = %5g   l2 = %5g  "
                    "F0 = %5g   F1 = %5g   F2 = %5g   going right\n",
                    l0, l1, l2, F0, F1, F2);

                l0=l1;
                l1=l2;
                F0=F1;
                F1=F2;
                l2=sg_lambda=fmin(1.25*l2, 0.95*lmax);
                if(!soah_shape_step_begin())
                    error_msg("%s: %4d: this should not happen.\n", __func__, __LINE__);
                copy_vec(init_u_c, init_u);
                copy_scal(init_p_c, init_p);
                F2=soah_shape_cost(init_u, init_p);
                continue;
            }
        }

        char message[1024];
        if(lmin_reached)
        {
            if(!soah_shape_step_acceptance_check(F0, Fs))
            {
                sprintf(message, "lambda = %e is too small.\n", l0);
                goto done_error;
            }
            else
            {
                sg_lambda=l0;
                Fs=F0;
                goto accept_step;
            }
        }

        if(lmin_reached)
        {
            if(!soah_shape_step_acceptance_check(F2, Fs))
            {
                sprintf(message, "lambda_max %g reached.\n", l2);
                goto done_error;
            }
            else
            {
                sg_lambda=l2;
                Fs=F2;
                goto accept_step;
            }
        }

        if(pit_found)
        {
            if(!soah_shape_step_acceptance_check(F1, Fs))
            {
                sprintf(message, "minimization failed.\n");
                goto done_error;
            }
            else
            {
                PRINTF("  pit:  l0 = %5g   l1 = %5g   l2 = %5g  "
                    "F0 = %5g   F1 = %5g   F2 = %5g   pit found\n",
                    l0, l1, l2, F0, F1, F2);
                sg_lambda=l1;
                Fs=F1;
                goto accept_step;
            }
        }

accept_step:
        {
            soah_shape_gradient(finl_v, finl_q);
            soah_shape_step_accept(init_u, init_p, finl_v, finl_q);
            soah_shape_step_termination_check();
            continue;
        }

done_error:
        {
            PRINTF("%s", message);

            AppendLogfile(saved_log);
            warning_msg("%s", message);
            DONE=TRUE;
            continue;
        }

    }
    soah_shape_optimize_end();
}

/****************************************************************************/
/*                                                                          */

/****************************************************************************/

void soah_shape_bloat(vector init_u, scalar init_p,
    vector finl_v, scalar finl_q)
{
    extern int _is_node_on_sigma(const grid3d_t *grid, int nn, int open_closed);
    extern double norm_sigma_scalar_zz(grid3d_t *grid, const scalar s1,
        const scalar s2);

    double h=grid->MAXH;

    PrintLog("%s: begin\n", __func__);

    soah_shape_optimize_begin(init_u, init_p);

    sg_lambda=0.0;
    if(!soah_shape_step_begin())
        exit(EXIT_FAILURE);
    clean_vec(init_u);
    clean_scal(init_p);
    clean_vec(finl_v);
    clean_scal(finl_q);
    PRINTF("Reference:\n");
    soah_shape_cost(init_u, init_p);
    soah_shape_gradient(finl_v, finl_q);

    double ref_cost=total_cost;
    double deriv=0.0;

    use_scalar(sg_xi)
    for(int i=0; i<grid->NN; i++)
    {
        if(_is_node_on_sigma(grid, i, FALSE))
        {
            deriv+=sg_xiF[i];
            sg_xiF[i]=1.0;
        }
        else
            sg_xiF[i]=0.0;
    }
    done_scalar(sg_xi)
    sg_norm=norm_sigma_scalar_zz(grid, sg_xi, NULL);
    PrintLog("       %10s   %-20s = %.8g\n", "deriv", "", deriv);


    shape_count++;

    sg_lambda=0.2*h;
    if(!soah_shape_step_begin())
        exit(EXIT_FAILURE);
    clean_vec(init_u);
    clean_scal(init_p);
    clean_vec(finl_v);
    clean_scal(finl_q);
    PRINTF("Deflated:\n");
    soah_shape_cost(init_u, init_p);
    AppendLogfile(saved_log);
    PrintLog("     Deflated:\n");
    PrintLog("       %10s   %-20s = %.8g\n", "total_cost", "", total_cost);
    double de_deriv=(total_cost-ref_cost)/sg_lambda;
    PrintLog("       %10s   %-20s = %.8g\n", "de_deriv", "", de_deriv);
    PRINTF("       de_deriv = %.8g\n", de_deriv);


    shape_count++;

    sg_lambda= -0.4*h;
    if(!soah_shape_step_begin())
        exit(EXIT_FAILURE);
    clean_vec(init_u);
    clean_scal(init_p);
    clean_vec(finl_v);
    clean_scal(finl_q);
    PRINTF("Inflated:\n");
    soah_shape_cost(init_u, init_p);
    AppendLogfile(saved_log);
    PrintLog("     Inflated:\n");
    PrintLog("       %10s   %-20s = %.8g\n", "total_cost", "", total_cost);
    double in_deriv=(total_cost-ref_cost)/sg_lambda;
    PrintLog("       %10s   %-20s = %.8g\n", "in_deriv", "", in_deriv);
    PRINTF("       in_deriv = %.8g\n", in_deriv);

    soah_shape_optimize_end();



}

/****************************************************************************
 *                                                                          *
 ****************************************************************************/

#include <time.h>

void soah_shape_optimize_test_gradient_at_point(vector init_u, scalar init_p,
    vector finl_v, scalar finl_q)
{
    double Fs, Gs_adj, Gs_lim;

    int test_node;

    PrintLog("%s: begin\n", __func__);

    soah_shape_optimize_begin(init_u, init_p);

    soah_shape_step_begin();
    clean_vec(init_u);
    clean_scal(init_p);
    clean_vec(finl_v);
    clean_scal(finl_q);
    Fs=soah_shape_cost(init_u, init_p);
    soah_shape_gradient(finl_v, finl_q);

    copy_vec(init_u, init_u_c);
    copy_scal(init_p, init_p_c);

    srand(clock());

    test_node= -1;
    for(int i=0; i<grid->NN; i++)
    {
        if(_is_node_on_sigma(grid, i, sigma_open)&&(!(rand()%17)))
        {
            test_node=i;
            break;
        }
    }

    if(test_node<0)
        error_msg("%s: test_node is negative\n", __func__);
    node_t *tn=grid->nodes+test_node;
    PrintLog("test_node = %d (%g, %g, %g)\n", test_node,
        tn->coords[0], tn->coords[1], tn->coords[2]);

    use_scalar(sg_xi);
    Gs_adj=sg_xiF[test_node];
    done_scalar(sg_xi);

    /*
    extern scalar sg_xi_weak;
    double Gs_adj_unfixed;
    use_scalar(sg_xi_weak)
    Gs_adj_unfixed=sg_xi_weakF[test_node];
    done_scalar(sg_xi_weak)
     */

    clean_scal(sg_xi);
    set_value_scal(sg_xi, test_node, -1.0);
    assemble_scal(sg_xi);

    set_reference_nodes(shape_nodes);
    // set_ref_nodes(alegrid, shape_nodes);
    AppendLogfile(saved_log);

    double lambda=grid->MAXH*0.5;
    sg_norm=1.0;

    for(int i=1; i<20; i++)
    {

        shape_count++;
        sg_lambda=lambda/(double)i;
        soah_shape_step_begin();
        clean_vec(init_u);
        clean_scal(init_p);
        double F2=soah_shape_cost(init_u, init_p);
        set_reference_nodes(shape_nodes);
        // set_ref_nodes(alegrid, shape_nodes);
        AppendLogfile(saved_log);
        Gs_lim=(F2-Fs)/sg_lambda;
        PrintLog("TEST: i = %d h = %g, Gs_adj = %g, Gs_lim = %g\n",
            i, sg_lambda, Gs_adj, Gs_lim);
        PrintLog("       %10s   %-20s = %.8g\n", "total_cost", "", total_cost);
        PRINTF("\rTEST: i = %d h = %g, Gs_adj = %g, Gs_lim = %g\n",
            i, sg_lambda, Gs_adj, Gs_lim);
    }


    soah_shape_optimize_end();

    PrintLog("%s: done\n", __func__);

}

